Pop nat analysis

1 DATA

#############################
#PERFORMANCE 
#############################

data_PERF <- import_data(dataset = "DATACOMPLET_PERF.csv", 
                         trait = "performance",
                         remove_testenvt = c("Grape","GF"), 
                         remove_pop = "WT3", 
                         remove_rate = NA)

## [1] "Data where the number of adults was higher than the number of eggs were not removed (85 and 1 tubes from 541 and 602 tubes for the first and thrid generation respectively)."

head(data_PERF)

##    Generation Experiment Original_environment   Population  Date_P Date_C_O
## 6          G0 Plasticity           Blackberry Blackberry45 3/10/18  4/10/18
## 8          G0 Plasticity           Blackberry Blackberry45 3/10/18  4/10/18
## 10         G0 Plasticity           Blackberry Blackberry45 3/10/18  4/10/18
## 12         G0 Plasticity           Blackberry Blackberry45 3/10/18  4/10/18
## 13         G0 Plasticity           Blackberry Blackberry45 3/10/18  4/10/18
## 14         G0 Plasticity           Blackberry Blackberry45 3/10/18  4/10/18
##      Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O Nb_adults Obs_A
## 6  19/10/2018  L1     C1    1       Blackberry       0    LO         0    LO
## 8  19/10/2018  L1     C3    1           Cherry       3    LO         2    LO
## 10 19/10/2018  L1     C5    1       Strawberry       0    LO         0    LO
## 12 19/10/2018  L1     C7    1       Strawberry       0    LO         0    LO
## 13 19/10/2018  L1     C8    1       Blackberry       1    LO         2    LO
## 14 19/10/2018  L1     C9    1           Cherry       2    LO         1    LO
##    EggScore EggScoreFive EggScoreSmall SA IndicG0 IndicG2 SAIndicG0
## 6         1            1             1  0       1       0         0
## 8         1            1             1  1       1       0         1
## 10        1            1             1  1       1       0         1
## 12        1            1             1  1       1       0         1
## 13        1            1             1  0       1       0         0
## 14        1            1             1  1       1       0         1
##                fruit_hab             fruit_hab_ng     fruit_gen       hab_gen
## 6  Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0 Blackberry_G0
## 8      Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
## 10 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
## 12 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
## 13 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0 Blackberry_G0
## 14     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
##            pop_gen      Rate
## 6  Blackberry45_G0       NaN
## 8  Blackberry45_G0 0.6666667
## 10 Blackberry45_G0       NaN
## 12 Blackberry45_G0       NaN
## 13 Blackberry45_G0 2.0000000
## 14 Blackberry45_G0 0.5000000

tapply(data_PERF$Nb_adults,list(data_PERF$Original_environment,data_PERF$Generation),length)

##             G0  G2
## Blackberry 241 320
## Cherry     258 143
## Strawberry  42 139

#############################
#EMERGENCE RATE
#############################

data_PERF_Rate <- import_data(dataset = "DATACOMPLET_PERF.csv", 
                         trait = "performance",
                         remove_testenvt = c("Grape","GF"), 
                         remove_pop = "WT3", 
                         remove_rate = TRUE)

## [1] "For the first generation, 85 tubes where the number of adults was higher than initial number of eggs were removed from 541 tubes."
## [1] "For the third generation, 1 tubes where the number of adults was higher than initial number of eggs were removed from 602 tubes."

###########################
#PREFERENCE 
###########################
data_PREF <- import_data(dataset = "DATACOMPLET_PREF.csv", 
                         trait = "preference",
                         remove_testenvt = NA, 
                         remove_pop = "WT3", 
                         remove_rate = NA)

head(data_PREF)

##   Generation Experiment BoxID     Date_P Original_environment   Population Line
## 1         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 2         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 3         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 4         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 5         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    2
## 6         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    2
##   Column Test_environment Nb_eggs   Date_C_O Obs_O SA IndicG0 IndicG2 SAIndicG0
## 1      1        Cranberry       0 13/12/2018    CD  1       1       0         1
## 2      2              Fig       0 13/12/2018    CD  1       1       0         1
## 3      3        Raspberry       0 13/12/2018    CD  1       1       0         1
## 4      4         Rosehips       0 13/12/2018    CD  1       1       0         1
## 5      1             Kiwi       0 13/12/2018    CD  1       1       0         1
## 6      2       Strawberry       1 13/12/2018    CD  1       1       0         1
##               fruit_hab             fruit_hab_ng     fruit_gen       hab_gen
## 1  Blackberry_Cranberry  Blackberry_Cranberry_G0 Blackberry_G0  Cranberry_G0
## 2        Blackberry_Fig        Blackberry_Fig_G0 Blackberry_G0        Fig_G0
## 3  Blackberry_Raspberry  Blackberry_Raspberry_G0 Blackberry_G0  Raspberry_G0
## 4   Blackberry_Rosehips   Blackberry_Rosehips_G0 Blackberry_G0   Rosehips_G0
## 5       Blackberry_Kiwi       Blackberry_Kiwi_G0 Blackberry_G0       Kiwi_G0
## 6 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
##           pop_gen
## 1 Blackberry33_G0
## 2 Blackberry33_G0
## 3 Blackberry33_G0
## 4 Blackberry33_G0
## 5 Blackberry33_G0
## 6 Blackberry33_G0

tapply(data_PREF$Nb_eggs,list(data_PREF$Original_environment,data_PREF$Generation),length)

##              G0   G2
## Blackberry  696 1176
## Cherry     1200  624
## Strawberry  252  480

###########################
#PREFERENCE 3 fruits
###########################
levels_test<-levels(data_PREF$Test_environment)
levels_original<-levels(data_PREF$Original_environment)


data_PREF_three <- import_data(dataset = "DATACOMPLET_PREF.csv", 
                         trait = "preference",
                         remove_testenvt = usefun::outersect(levels_test, 
                                                             levels_original), 
                         remove_pop = "WT3", 
                         remove_rate = NA)

head(data_PREF_three)

##    Generation Experiment BoxID     Date_P Original_environment   Population
## 6          G0 Plasticity   860 19/09/2018           Blackberry Blackberry33
## 9          G0 Plasticity   860 19/09/2018           Blackberry Blackberry33
## 11         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33
## 15         G0 Plasticity   884 19/09/2018           Blackberry Blackberry33
## 16         G0 Plasticity   884 19/09/2018           Blackberry Blackberry33
## 20         G0 Plasticity   884 19/09/2018           Blackberry Blackberry33
##    Line Column Test_environment Nb_eggs   Date_C_O Obs_O SA IndicG0 IndicG2
## 6     2      2       Strawberry       1 13/12/2018    CD  1       1       0
## 9     3      1       Blackberry       0 13/12/2018    CD  0       1       0
## 11    3      3           Cherry       0 13/12/2018    CD  1       1       0
## 15    1      3           Cherry       0 13/12/2018    CD  1       1       0
## 16    1      4       Strawberry       0 13/12/2018    CD  1       1       0
## 20    2      4       Blackberry       0 13/12/2018    CD  0       1       0
##    SAIndicG0             fruit_hab             fruit_hab_ng     fruit_gen
## 6          1 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0
## 9          0 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0
## 11         1     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0
## 15         1     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0
## 16         1 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0
## 20         0 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0
##          hab_gen         pop_gen
## 6  Strawberry_G0 Blackberry33_G0
## 9  Blackberry_G0 Blackberry33_G0
## 11     Cherry_G0 Blackberry33_G0
## 15     Cherry_G0 Blackberry33_G0
## 16 Strawberry_G0 Blackberry33_G0
## 20 Blackberry_G0 Blackberry33_G0

tapply(data_PREF_three$Nb_eggs,list(data_PREF_three$Original_environment,
                                    data_PREF_three$Generation),length)

##             G0  G2
## Blackberry 174 294
## Cherry     300 156
## Strawberry  63 120

dim(data_PREF_three)

## [1] 1107   21

dim(data_PREF)

## [1] 4428   21

2 NUMBER OF EGGS

2.1 Exploration data

tapply(data_PERF$Nb_eggs, list(data_PERF$Original_environment, data_PERF$Test_environment, data_PERF$Generation), mean, na.rm = TRUE)

## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry   4.795181  7.831169   4.649351
## Cherry      21.453488 18.395349  10.694118
## Strawberry  19.785714 24.000000  14.428571
## 
## , , G2
## 
##            Blackberry   Cherry Strawberry
## Blackberry   103.9720 121.8491   98.71963
## Cherry       131.7174 133.8400  103.51064
## Strawberry   119.1739 119.4894  111.56522

ggplot2::ggplot(data = data_PERF[data_PERF$Generation=="G0",], 
                aes(x = Test_environment, y = Nb_eggs, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF, 
                aes(x = Nb_eggs, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

2.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))

## [1] 1.294236

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 

## [1] 0.3379054

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5, 1]))

## [1] 1.38122

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[5, 1]) )

## [1] 0.3246791

## F test for SA
(Fratio_NonGen <- (anova(m2)[4,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 1.243819

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[6, 1]) )

## [1] 0.3460299

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1]))))

## [1] 0.08701221

(r2_SA_nongenet <- 1-(anova(m2)[6, 3]/((anova(m2)[4, 2]+anova(m2)[6, 2])/(anova(m2)[4, 1]+anova(m2)[6, 1]))))

## [1] 0.05745281

##Plot

(PLOT_eggs_G0 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_eggs", gen = "G0"))

(PLOT_eggs_G2 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_eggs", gen = "G2"))

3 NUMBER OF ADULTS

3.1 Exploration data

tapply(data_PERF$Nb_adults, list(data_PERF$Original_environment, data_PERF$Test_environment, data_PERF$Generation), mean, na.rm = TRUE)

## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry   3.380952  3.443038   2.641026
## Cherry       6.383721  5.372093   2.197674
## Strawberry  11.357143 11.428571   4.642857
## 
## , , G2
## 
##            Blackberry   Cherry Strawberry
## Blackberry   27.31776 23.18868   13.15888
## Cherry       22.91304 23.32000   11.53191
## Strawberry   25.65217 22.36170   16.04348

ggplot2::ggplot(data = data_PERF[data_PERF$Generation=="G0",], 
                aes(x = Test_environment, y = Nb_adults, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF[data_PERF$Generation=="G2",], 
                aes(x = Test_environment, y = Nb_adults, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF, 
                aes(x = Nb_adults, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

3.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))

## [1] 2.08912

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 

## [1] 0.2441119

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5, 1]))

## [1] 2.407889

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[5, 1]) )

## [1] 0.2185257

## F test for SA
(Fratio_NonGen <- (anova(m2)[4,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 1.152247

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[6, 1]) )

## [1] 0.3617404

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1]))))

## [1] 0.2603399

(r2_SA_nongenet <- 1-(anova(m2)[6, 3]/((anova(m2)[4, 2]+anova(m2)[6, 2])/(anova(m2)[4, 1]+anova(m2)[6, 1]))))

## [1] 0.03666627

3.3 Analysis with eggs as covariable

#######################################################
## Should we consider the number of eggs?   ###
#######################################################

# Original
m1 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF)


## Correction for number of eggs
m2 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            log(Nb_eggs+1), 
          data = data_PERF)


## With egg score
m3 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScore, 
          data = data_PERF)

## Compare with 5 egg scores
m4 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreFive, 
          data = data_PERF)

## Compare with EggScoreSmall
m5 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreSmall, 
          data = data_PERF)



MuMIn::model.sel(m1, m2, m3, m4, m5)

## Model selection table 
##     (Int) hab_gen pop_gen SA IG0:SA Org_env:Tst_env IG0:Org_env:Tst_env
## m2 0.3032       +       +  +      +               +                   +
## m5 0.7887       +       +  +      +               +                   +
## m3 0.7890       +       +  +      +               +                   +
## m4 0.7759       +       +  +      +               +                   +
## m1 1.0270       +       +  +      +               +                   +
##    log(Nb_egg+1) EgS ESF ESS             family df    logLik   AICc  delta
## m2         0.472             gaussian(identity) 63 -1064.139 2261.8   0.00
## m5                         + gaussian(identity) 67 -1155.928 2454.4 192.59
## m3                 +         gaussian(identity) 65 -1161.139 2460.3 198.50
## m4                     +     gaussian(identity) 66 -1160.686 2461.6 199.84
## m1                           gaussian(identity) 62 -1248.932 2629.1 367.31
##    weight
## m2      1
## m5      0
## m3      0
## m4      0
## m1      0
## Models ranked by AICc(x)

# Models are not all fitted to the same data: because 6 tubes without Nb_eggs are missing for m2, m3, m4 and m5


## Cl= The best model is when the number of eggs is considered as a continuous variable
### model m2 provides a better description of the data than model m1 






#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment + log(Nb_eggs+1),
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5, 1]))

## [1] 9.32529

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[5, 1])) 

## [1] 0.0552651

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 + 
            log(Nb_eggs+1), 
          data = data_PERF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 12.31994

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )

## [1] 0.03919852

## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))

## [1] 0.876109

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )

## [1] 0.4183184

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))

## [1] 0.7389023

(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))

## [1] -0.03196271

##Plot

(PLOT_adult_G0 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_adults", gen = "G0"))

(PLOT_adult_G2 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_adults", gen = "G2"))

4 EMERGENCE RATE

4.1 Exploration data

tapply(data_PERF_Rate$Rate, list(data_PERF_Rate$Original_environment, 
                                 data_PERF_Rate$Test_environment, 
                                 data_PERF_Rate$Generation), mean, na.rm = TRUE)

## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry  0.6125868 0.4703813  0.5248347
## Cherry      0.2420863 0.3199467  0.2176492
## Strawberry  0.2888558 0.3595403  0.3373855
## 
## , , G2
## 
##            Blackberry    Cherry Strawberry
## Blackberry  0.2579090 0.1923140  0.1345134
## Cherry      0.1855892 0.1801664  0.1103531
## Strawberry  0.2206513 0.1950778  0.1619433

ggplot2::ggplot(data = data_PERF_Rate[data_PERF_Rate$Generation=="G0",], 
                aes(x = Test_environment, y = Rate, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF_Rate[data_PERF_Rate$Generation=="G2",], 
                aes(x = Test_environment, y = Rate, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF_Rate, 
                aes(x = Rate, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

lattice::xyplot(Rate~Nb_eggs|Original_environment*Test_environment, 
                data=data_PERF_Rate)

lattice::xyplot(Rate~EggScore|Original_environment*Test_environment, 
                data=data_PERF_Rate)

lattice::bwplot(Rate~EggScore|Original_environment*Test_environment, 
                data=data_PERF_Rate)

lattice::bwplot(Rate~EggScoreFive|Original_environment*Test_environment, 
       data=data_PERF_Rate)

lattice::bwplot(Rate~EggScoreSmall|Original_environment*Test_environment, 
       data=data_PERF_Rate)

AvEmergenceRate <- tapply(data_PERF_Rate$Rate, 
                          list(data_PERF_Rate$EggScoreFive,
                              data_PERF_Rate$Original_environment,
                              data_PERF_Rate$Test_environment),mean)
tapply(data_PERF_Rate$Rate, list(data_PERF_Rate$EggScoreFive,
                                 data_PERF_Rate$Original_environment,
                                 data_PERF_Rate$Test_environment), length)

## , , Blackberry
## 
##   Blackberry Cherry Strawberry
## 1         74     73         14
## 2         42     22         15
## 3         41     17         20
## 4         10     10          8
## 5          3      4          2
## 
## , , Cherry
## 
##   Blackberry Cherry Strawberry
## 1         59     72         12
## 2         31     18         19
## 3         50     21         21
## 4         23     11          9
## 5          2      7         NA
## 
## , , Strawberry
## 
##   Blackberry Cherry Strawberry
## 1         60     80         20
## 2         55     24         13
## 3         39     16         16
## 4          7      6          9
## 5         NA     NA          2

AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Blackberry"][, "Blackberry"]

##         1         2         3         4         5 
##       NaN 1.0258692 0.7334304 0.9735309 0.1709919

AvEmergenceRate[, , "Cherry"][, "Blackberry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]

##         1         2         3         4         5 
##       NaN 0.8562518 1.0861981 0.8514873 1.4918117

AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]

##         1         2         3         4         5 
##       NaN 1.0522488 1.0774054 0.9037047 0.4651315

4.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF_Rate)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))

## [1] 30.89092

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 

## [1] 0.0114893

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF_Rate)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5, 1]))

## [1] 59.30375

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[5, 1]) )

## [1] 0.004550853

## F test for SA
(Fratio_NonGen <- (anova(m2)[4,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 3.305405

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[6, 1]) )

## [1] 0.1666402

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1]))))

## [1] 0.9357984

(r2_SA_nongenet <- 1-(anova(m2)[6, 3]/((anova(m2)[4, 2]+anova(m2)[6, 2])/(anova(m2)[4, 1]+anova(m2)[6, 1]))))

## [1] 0.3656236

4.3 Analysis with eggs as covariable

#######################################################
## Should we consider the number of eggs?   ###
#######################################################

# Original
m1 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF_Rate)


## Correction for number of eggs
m2 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            log(Nb_eggs+1), 
          data = data_PERF_Rate)


## With egg score
m3 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScore, 
          data = data_PERF_Rate)

## Compare with 5 egg scores
m4 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreFive, 
          data = data_PERF_Rate)

## Compare with EggScoreSmall
m5 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreSmall, 
          data = data_PERF_Rate)



MuMIn::model.sel(m1, m2, m3, m4, m5)

## Model selection table 
##     (Int) hab_gen pop_gen SA IG0:SA Org_env:Tst_env IG0:Org_env:Tst_env
## m2 0.7356       +       +  +      +               +                   +
## m1 0.6018       +       +  +      +               +                   +
## m3 0.6200       +       +  +      +               +                   +
## m4 0.6210       +       +  +      +               +                   +
## m5 0.6186       +       +  +      +               +                   +
##    log(Nb_egg+1) EgS ESF ESS             family df   logLik  AICc delta weight
## m2      -0.07341             gaussian(identity) 63 -134.356 403.4  0.00      1
## m1                           gaussian(identity) 62 -149.207 430.9 27.42      0
## m3                 +         gaussian(identity) 65 -146.871 433.0 29.61      0
## m4                     +     gaussian(identity) 66 -146.858 435.3 31.88      0
## m5                         + gaussian(identity) 67 -146.342 436.6 33.15      0
## Models ranked by AICc(x)

## Cl= The best model is when the number of eggs is considered as a continuous variable
### model m2 provides a better description of the data than model m1 







#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment +
            log(Nb_eggs+1),
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF_Rate)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5, 1]))

## [1] 25.99987

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[5, 1])) 

## [1] 0.01458567

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            log(Nb_eggs+1), 
          data = data_PERF_Rate)


## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 34.75823

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )

## [1] 0.009741755

## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))

## [1] 2.109581

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )

## [1] 0.2423154

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))

## [1] 0.8940628

(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))

## [1] 0.217157

##Plot

(PLOT_rate_G0 <- plot_RTP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G0"))

(PLOT_rate_G2 <- plot_RTP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G2"))

5 PREFERENCE

5.1 Twelve fruits

5.1.1 Exploration data

tapply(data_PREF$Nb_eggs, list(data_PREF$Original_environment, 
                                 data_PREF$Test_environment, 
                                 data_PREF$Generation), mean, na.rm = TRUE)

## , , G0
## 
##              Apricot Blackberry Blackcurrant    Cherry Cranberry       Fig
## Blackberry 0.2758621  0.7931034    0.1379310 0.4137931 0.2068966 0.2586207
## Cherry     0.3300000  0.5000000    1.2600000 2.4800000 0.3800000 0.7400000
## Strawberry 0.5714286  1.1904762    0.6666667 1.6666667 0.5714286 0.4285714
##                Grape      Kiwi Raspberry  Rosehips Strawberry    Tomato
## Blackberry 0.2068966 0.0862069 0.4655172 0.3793103  0.5862069 0.2241379
## Cherry     0.9700000 1.0500000 1.5800000 1.2400000  0.6700000 0.6969697
## Strawberry 0.4761905 0.1428571 0.6666667 0.3809524  0.8571429 0.0952381
## 
## , , G2
## 
##              Apricot Blackberry Blackcurrant   Cherry Cranberry      Fig
## Blackberry  6.826531   24.62245     18.61224 23.41837  8.428571 12.73469
## Cherry     10.423077   16.09615     19.09615 30.38462  7.480769 13.46154
## Strawberry 15.450000   21.92500     22.12500 24.92500  7.200000 13.02500
##               Grape     Kiwi Raspberry Rosehips Strawberry   Tomato
## Blackberry 24.55102 16.57143  16.59184 11.23469   6.938776 12.53061
## Cherry     11.67308 11.90385  12.78846 11.46154  10.673077 10.59615
## Strawberry 12.55000 14.87500  19.45000 16.20000  12.775000 14.90000

ggplot2::ggplot(data = data_PREF, 
                aes(x = Test_environment, y = Nb_eggs, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PREF, 
                aes(x = Nb_eggs, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

5.1.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PREF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))

## [1] 10.1742

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 

## [1] 0.004407544

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            BoxID, 
          data = data_PREF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 9.833047

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )

## [1] 0.004993505

## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))

## [1] 2.500166

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )

## [1] 0.1287796

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))

## [1] 0.2864799

(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))

## [1] 0.0638364

5.2 Three fruits

5.2.1 Exploration data

tapply(data_PREF_three$Nb_eggs, list(data_PREF_three$Original_environment, 
                                 data_PREF_three$Test_environment, 
                                 data_PREF_three$Generation), mean, na.rm = TRUE)

## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry  0.7931034 0.4137931  0.5862069
## Cherry      0.5000000 2.4800000  0.6700000
## Strawberry  1.1904762 1.6666667  0.8571429
## 
## , , G2
## 
##            Blackberry   Cherry Strawberry
## Blackberry   24.62245 23.41837   6.938776
## Cherry       16.09615 30.38462  10.673077
## Strawberry   21.92500 24.92500  12.775000

ggplot2::ggplot(data = data_PREF_three, 
                aes(x = Test_environment, y = Nb_eggs, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PREF_three, 
                aes(x = Nb_eggs, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

5.2.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment +
            BoxID,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PREF_three)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5, 1]))

## [1] 18.38385

(pvalue <- 1 - pf(Fratio, 1, anova(m1)[5, 1])) 

## [1] 0.02331823

#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            BoxID, 
          data = data_PREF_three)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))

## [1] 12.77418

(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )

## [1] 0.0374429

## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))

## [1] 3.536493

(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )

## [1] 0.1566089

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))

## [1] 0.746421

(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))

## [1] 0.3880511

5.2.3 Plot

PLOT_pref_G0 <- plot_RTP_residuals(dataset = data_PREF_three, trait = "Nb_eggs", gen = "G0")
PLOT_pref_G2 <- plot_RTP_residuals(dataset = data_PREF_three, trait = "Nb_eggs", gen = "G2")

6 PLOT

legend <- lemon::g_legend(PLOT_pref_G0)

############## FIRST GENERATION 
LOCAL_ADAPTATION_PLOT <- cowplot::ggdraw() +
  cowplot::draw_plot(PLOT_eggs_G0+theme(legend.position = "none", 
                                        plot.title = element_blank()),
            x =0, y = 0.5, width = 0.36, height = 0.46) +
  cowplot::draw_plot(PLOT_adult_G0+theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.4, y = 0.5, width =0.36, height = 0.46) +
  cowplot::draw_plot(PLOT_rate_G0+theme(legend.position = "none", 
                                        plot.title = element_blank()),
            x =0, y = 0, width = 0.36, height = 0.46) +
  cowplot::draw_plot(PLOT_pref_G0+theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.4, y = 0, width = 0.36, height = 0.46) +
  cowplot::draw_plot(legend, x = 0.85, y = 0.5, width = 0.0001, height = 0.0001) +
  cowplot::draw_plot_label(c("A", "B", "C", "D"), 
                           c(0.01, 0.4,0.01, 0.4), 
                  c(1, 1,0.5,0.5), size = 19)
# + cowplot::draw_plot_label("First generation", x = 0.2, y = 1 , size = 14) 
LOCAL_ADAPTATION_PLOT

cowplot::save_plot(file =here::here("figures","Reciprocal_experiment_FirstGeneration.pdf"),
                  LOCAL_ADAPTATION_PLOT, base_height = 20/cm(1),
                  base_width = 30/cm(1), dpi = 1200)


############## THIRD GENERATION 
LOCAL_ADAPTATION_PLOT_THIRD <- cowplot::ggdraw() +
  cowplot::draw_plot(PLOT_eggs_G2+theme(legend.position = "none", 
                                        plot.title = element_blank()),
            x =0, y = 0.5, width = 0.36, height = 0.46) +
  cowplot::draw_plot(PLOT_adult_G2+theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.4, y = 0.5, width =0.36, height = 0.46) +
  cowplot::draw_plot(PLOT_rate_G2+theme(legend.position = "none", 
                                        plot.title = element_blank()),
            x =0, y = 0, width = 0.36, height = 0.46) +
  cowplot::draw_plot(PLOT_pref_G2+theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.4, y = 0, width = 0.36, height = 0.46) +
  cowplot::draw_plot(legend, x = 0.85, y = 0.5, width = 0.0001, height = 0.0001) +
  cowplot::draw_plot_label(c("A", "B", "C", "D"), 
                           c(0.01, 0.4,0.01, 0.4), 
                  c(1, 1,0.5,0.5), size = 19)
# + cowplot::draw_plot_label("First generation", x = 0.2, y = 1 , size = 14) 
LOCAL_ADAPTATION_PLOT_THIRD

cowplot::save_plot(file =here::here("figures","Reciprocal_experiment_ThirdGeneration.pdf"),
                  LOCAL_ADAPTATION_PLOT_THIRD, base_height = 20/cm(1),
                  base_width = 30/cm(1), dpi = 1200)